Growing length and time scales in a suspension of athermal particles 
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We simulate a relaxation process of noii-brownian particles in a sheared viscous medium; the 
small shear strain is initially applied to a system, which then undergoes relaxation. The relaxation 
time and the correlation length are estimated as functions of density, which algebraically diverge at 
the jamming density. This implies that the relaxation time can be scaled by the correlation length 
using the dynamic critical exponent, which is estimated as 4.6(2). It is also found that shear stress 
undergoes power-law decay at the jamming density, which is reminiscent of critical slowing down. 
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In general (thermal or athermal) particulate systems at 
high density, structural rearrangement of the constituent 
particles is difficult due to the exclusion volume effect 
so that the structural relaxation time and the viscosity 
drastically increase. In particular, zero temperature sys- 
tems such as granular materials or emulsions acquire the 
elastic moduli above a certain density [l|, \^ . This rigidity 
transition, which is also referred to as the jamming tran- 
sition, is accompanied by some power-law behaviors that 
are characteristic to critical phenomena: e.g. the grow- 
ing correlation leiigth in terms of spatially heterogeneous 
diffusion [2, ^, la, lo, LZf . Note that, indeed prior to those 
in athermal systems, such dynamical heterogeneities are 
widely observed in thermal systems such as supercooled 
liquids and dense colloids [8J, |9|, [lOl, [ill . This suggests the 
potential relation between the glass a nd j amming transi- 
tions, although still controversial [12Lll3l.ll4ll5l|. 

Aside from the connection with glass transitions, how- 
ever, the nature of the jamming transition itself is still not 
clear. Provided that jamming is a critical phenomenon, 
the critical exponents that describe the divergence of the 
length and the time scales play an essential role in clarify- 
ing the underlying mechanism that yields the dynamical 
heterogeneity and classifying the jamming transition into 
a universality class (if any). There are several attempts 
to estimate such exponents, in particular that for the 
correlation length; £, ^ \4> — (t)j\~^ , where </> denotes the 
density and (j)j is the critical density. Finite-size scaling 
reveals that z/ ~ 0.7 for both two and three dimensional 
systems, although the correlation length is not explicitly 
analyzed [2]. Later, the correlation length is defined us- 
ing the concept of the dynamical heterogeneity so that 
the exponent v is estimated [5|, la, lZ|, [la, ll7|. Note that, 
however, the above attempts involve two dimensional sys- 
tems and there is no direct estimation of the exponent v 
in three dimensional systems. Recalling that the dimen- 
sionality plays an essential role in conventional critical 
phenomena, one must investigate the three dimensional 
case. 



Another important quantity that characterizes a crit- 
ical phenomenon is the characteristic time. Wyart et 
al. derive the characteristic frequency of jammed sys- 



tems using the normal mode analysis. They obtain 
LOc '-^ {4> — 4'j)'' : where C = 0.5 [l8|. Interestingly, their 
discussion does not depend on the dimensionality of the 
system. However, the characteristic time in unjammed 
systems is still not estimated. Along the line of thought, 
in this paper, we investigate the behaviors of the relax- 
ation time and the correlation length in a three dimen- 
sional unjammed system at zero temperature. It is found 
that the relaxation time r and the correlation length ^ 
increases obeying power laws with respect to the density; 
T ^ (</>,/ — 0)~'^'^ and^ ~ {4>j — (t))~^-'^. These results lead 
to z ~ 4.6, which coincides with that of a Lennard- Jones 
glass [H. 

We consider macroscopic particles in a viscous medium 
so that temperature does not play any role. We neglect 
hydrodynamic and electric interactions because we are 
interested in the nature of the jamming critical point. 
The particles arc monodisperse, the diameter of which is 
denoted by d. Note that the particles are elastic so that 
the interaction between particles i and j is described by 
linear repulsive force; i.e., fij 



khijUij, where k denotes 



the elastic constant, n.; 



Cij\, and hij 



d- |r,j| 
if d < 



denoting the overlap length. (Note that hij 
I'Tijl)- We consider only unjammed systems, the volume 
density of which is less than the jamming density; </> < 
(f>j ~ 0.639 J4|. Unless otherwise indicated, the system 
contains 6.4 x 10^ particles. 

A relaxation process is realized in such a way that an 
equilibrated system is perturbed at i = and then un- 
dergoes time evolution. To prepare an initial equilibrated 
system under periodic boundary conditions, we adopt 
the conjugate gradient method, by which randomly- 
distributed particles relax to a zero energy state (no 
overlaps between particles). Then we perturb this equi- 
librated system by applying pure shear with respect to 
the {y,z) plane; i.e., each particle is instantaneously dis- 
placed by the following Affine deformation. 



Vi + ezi 
Zi + e.yi 



(1) 
(2) 



where e denotes the shear strain. At the same time, pe- 
riodic boundary conditions are slightly modified in order 




FIG. 1: (a) The relaxation of shear stress at various densities 
(see the legends) . Lower curves correspond to lower densities. 
The relaxation time grows as the density increases towards 
the critical density, (b) Rescaled shear stress S/S as a func- 
tion of rescaled time, t/T{(t>), where t{<I)) is the characteristic 
time that depends on the density, and S is approximately 



to adapt the shear strain; the adjacent imaginary cells 
with respect to the y and z directions are displaced by 
eL [20|, while an ordinary periodic boundary condition 
is used along the x direction. Due to the shear strain, 
overlaps between particles appear so that the system ac- 
quires nonzero elastic energy at t = 0. Then the sys- 
tem begins to relax and eventually reaches a new stable 
state of zero energy. In order to mimic the dynamics 
of macroscopic particles in a highly viscous medium, we 
adopt overdamp dynamics; i.e., jVi = ^^(ij- Through- 
out this study, we adopt the units in which the mass, 
the diameter, and the mobility 7 are unity. The elastic 
constant k is set to be 7^/to. This procedure may be 
realizable in experiments using macroscopic particles in 
a viscous medium 21[ , where the shear strain is applied 
via the viscous medium and the gravity can be canceled 
by density-matching. 

First, we discuss the relaxation of macroscopic quan- 
tities, in particular shear stress, which is defined via 
the virial [22]. Because the temporal behavior of other 
macroscopic quantities such as pressure and energy are 
essentially the same, we focus shear stress here. The re- 
laxation of shear stress at each density is shown in FIG. 
[T] (a), where the initial strain e is 1.0 x 10^'^. We also 
test e — 0.01 and e — 0.05 to find that the result does 
not qualitatively depends on the initial strain. It is im- 
portant to remark that several samples at each density 
show quantitatively the same relaxation behavior; less 
than ±10 % in the characteristic time defined below. 

Then we define the characteristic time from the re- 
laxation behaviors shown in FIG. [T] (a). These tempo- 
ral behaviors indeed collapse by rescaling the time with 







FIG. 2: The relaxation time as a function of density. Solid line 
in each panel is proportional to r ~ $~ ' , where <E> = 0,/ — </!> 
and (j>j = 0.639. (a) The relaxation time normalized by tq: 
To = 2 for e = 0.05 , To = 1 otherwise. The power-law diver- 
gence and the exponent do not depend on the initial strain, 
(b) There is no system-size dependence of the relaxation time. 
The initial strain e = 0.05. 



the characteristic time at each density, t(0), as shown 
in Fig. [T] (b). Note that shear stress is also rescaled by 
S ~ 0.5So/\/t, which indicates that the temporal behav- 
ior of shear stress is described by 



S{t) ^ Sot 



-»g-t/T 



(3) 



with a ~ 0.5. By fitting the shear stress relaxation with 
Eq. Q, we obtain a = 0.55(5). The relaxation time 
t(0) drastically increases at higher densities as shown in 
FIG. [H (a), which indicates the power-law divergence of 
the relaxation time at the jamming density. 



T0( 



-C 



(4) 



where tq denotes the time constant that does not de- 
pend on the density and the exponent is estimated as 
C = 3.3(1). However, note that this exponent is valid 
only in higher density region, </> > 0.60. There seems to 
be a crossover to the different behavior in lower density 
region, </> < 0.60, where ( = 1.5(1) if we assume Eq. ^. 
Note that this is not due to a finite-size effect as illus- 
trated in FIG. [H (b), which indicates the relaxation time 
for three different systems: TV = 1000, 8000, 64000. This 
is indeed due to the qualitative change in the particle 
dynamics, because the correlation function discussed be- 
low, Eq. ([7]), cannot detect any cooperative motion in 
the lower density region, < 0.60. However, as we focus 
on the critical nature of the jamming transition here, we 
do not further discuss this crossover. 

From Eqs. ([T|) and (g]), it is expected that the relax- 
ation of shear stress is described by a simple power law at 
the jamming density, S{t) ^ Sot~'^. Figure [3] (a) shows 
the stress relaxation at </> = 0.637, which obeys a power- 
law decay for a considerable duration (up to six orders 
of magnitude). We remark that this power-law relax- 
ation, which is quite similar to critical slowing down, is 
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FIG. 3: (a) The relaxation of shear stress for (j> = 0.637 and 
e = 0.05. The solid hne is proportional to t""'^. (b) The 
relaxation of the average coordination number for (p — 0.637 
and e = 0.05. (c) The relaxation of the average magnitude 
of interparticle force, {/). Time is rescaled by t((/>) shown in 
FIG. [3 (a). Note that {/) is also rescaled by / ~ r°^ The 
initial strain is 0.001. 
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FIG. 4: (a) The collapse of the correlation function defined by 
Eq. ([Xl . Note that the length is normalized by the correlation 
length, ^. (b) The correlation length as a function of <E> = 
(l>.j — (p, where (pj = 0.639. The dashed line is proportional 
to $^''-^^. Note that the correlation length is normalized by 
^0, which depends on the initial strain e. Here ^o is set to be 
unity for ^o = 0.05, while ^o = 0.5 for e = 0.001 and ^o = 0.8 
for e = 0.01. 



consistent with a theory in which shear stress is an or- 
der parameter that undergoes marginal stability at the 
jamming transition point [23l. |24|. 

The slow relaxation is also apparent in view of the 
particle dynamics. A quantity that involves the stresses 
and energy is the magnitude of interparticle force. We 
thus define the average magnitude of interparticle force 
as 



(/) 



2 



Eif'^i' 



(5) 



i>j 



where Z is the coordination number. We remark that Z is 
almost constant during the relaxation process as shown in 
FIG. [3] (b), while the average magnitude of force relaxes 
as 



(/) « r 



-/3p-*/r 



(6) 



where /3 ~ 0.3. The relaxation profile of (/) is shown in 
FIG. [3] (c), where the time is rescaled by the relaxation 
time T defined by Eq. ([3]). Furthermore, we observe 
that energy and pressure relax in essentially the same 
manner; i.e., P oc t~^e~*^^ and E oc t~'^^e~^*^^. This 
is indeed trivial because P - 2{NZ)-^{f)d/V and E - 
2{N Z)^^ {P) / kV , where d and k denote the diameter 
and the stiffness of the particles, respectively. We also 
remark that / obeys an exponential distribution so that 

iP) - if)'- 

We then investigate the spatial correlation of the par- 
ticle motion during the relaxation process. To this end, 
we compare the initial and the final configurations to 
define the displacement vector of each particle; Ax,; = 



,(1) 



,(") 



A(x, 



(0)^ 



where xj ' and x,- ' 



are the initial 



and the final positions of particle i, and A(r) is the Affine 
deformation vector at position r, represented by Eqs. ^ 
and (21). Using the displacement vectors Ax^, we define 
the following correlation function [9|, |25| . 



G(r) 



E.>,Ax,.Ax,<5(r-|xr^-x^^|) 



(0) 



,(0)i 



E.>,^(^-|x: 



(0) 



M, 



(7) 



This quantifies the extent of chainlike cooperative mo- 
tion, which is universally observed in dense particulate 
systems; i.e., supercooled liquids, colloids, and grains 
[i, m, li^. Hi]. Figure H (a) shows the collapse of the 



correlation function, where ^G{r/£,) is plotted as a func- 
tion of r/^. This indicates that ^ is the correlation length 
and the correlation function G{r) is approximately expo- 
nential. As is shown in FIG.|l](b), this correlation length 
f increases as the density approaches the jamming den- 
sity; ^ ^ S,o{4>J — 4')~'' ^ where v ~ 0.7. This value is 
consistent with the shift exponent for the static rigid- 
ity transition in three dimensions [2]. Also, within the 
present numerical accuracy, it is indistinguishable from 
the exponent obtained from the dynamical heterogeneity 
in two dimensional systems [SJ, y, [sj, |6|, lD, 116| 

From the present numerical results we can estimate the 
dynamic critical exponent, z = (,lv. We replot the cor- 
relation length and the relaxation time in FIG. [51 where 
we can estimate z as 4.6(2). This value is quite differ- 
ent from that obtained in two-dimensional air-fluidized 
beads, where z ~ 1.4 y. Although the difference may 
be due to the dimensionality, we do not have any definite 
explanation on this difference at this point. The struc- 
tural relaxation process may strongly depend on the spe- 
cific dynamics, particularly an agitation method by which 
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FIG. 5: The correlation length (normalized by ^o) dependence 
of the relaxation time (normalized by to). Data are taken 
from FIGS. [2] and m The dashed line is proportional to ^'^■''. 



energy is injected to the system. Nevertheless, it is note- 
worthy that almost the same dynamic critical exponent 
{z ~ 4.55 ± 0.2) is found in a binary Lennard-Jones (LJ) 
supercooled liquid, where the correlation length and the 
characteristic time for the structural relaxation diverges 
towards zero temperature [19[. This result, together with 
a recent study [15|, implies that a binary LJ glass is dom- 
inated by a zero temperature critical point, which may 
have close relation to the jamming transition. Although 
the correspondence between an LJ glass and the present 
system is not apparent, this coincidence may provide an- 
other evidence for the close relation between the glass 
and jamming transitions. 

To summarize, we show power-law divergence of the 
relaxation time and the correlation length as the den- 
sity approaches the jamming density from below; t/tq ^ 
(0j - cl))-'^ ~ (C/Co)^ where C = 3.3(1) and z = 4.6(2). 
The correlation length is defined in terms of stringlike 
cooperative motion, which is widely observed in dense 
particulate systems. It is also found that the stress re- 
laxation near the jamming density obeys critical slowing 
down. Although the result is restricted to a three dimen- 
sional unjammed system, we remark that the length and 
time scales in jammed systems (where (f) > 4>,j) and the 
effect of dimensionality is currently investigated by using 
the present relaxation method. 
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